{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "9a538b6f-df3c-4d30-b89f-cafb5e2597e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "24366129-a37a-4007-8d33-121d4385a932",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('/zfs/disinfo/dcweekly_contrast/story_text/dc_weekly_topic_scores.csv')\n",
    "# Make a binary indicator for whether we are in the AI-generated text period of time\n",
    "df['period'] = df.month.apply(lambda x: 1- int(x=='2023-06'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "becd539f-5d61-4c03-b6b4-13efe4963227",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                   guns   R-squared:                       0.374\n",
      "Model:                            OLS   Adj. R-squared:                  0.374\n",
      "Method:                 Least Squares   F-statistic:                     2017.\n",
      "Date:                Tue, 23 Apr 2024   Prob (F-statistic):               0.00\n",
      "Time:                        08:05:56   Log-Likelihood:                -470.34\n",
      "No. Observations:               10144   AIC:                             948.7\n",
      "Df Residuals:                   10140   BIC:                             977.6\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "period         0.1946      0.007     26.513      0.000       0.180       0.209\n",
      "Israel         0.2840      0.007     39.771      0.000       0.270       0.298\n",
      "Ukraine        0.3399      0.007     50.799      0.000       0.327       0.353\n",
      "intercept      0.1745      0.006     27.118      0.000       0.162       0.187\n",
      "==============================================================================\n",
      "Omnibus:                       41.450   Durbin-Watson:                   1.869\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               41.915\n",
      "Skew:                          -0.156   Prob(JB):                     7.91e-10\n",
      "Kurtosis:                       2.962   Cond. No.                         5.35\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# Perform linear regression\n",
    "# Add a constant term for the intercept\n",
    "df['intercept'] = 1\n",
    "\n",
    "# Specify the model\n",
    "independent_vars = ['period', 'Israel', 'Ukraine', 'intercept']\n",
    "dependent_var = 'guns'\n",
    "\n",
    "# Fit the model\n",
    "guns_model = sm.OLS(df[dependent_var], df[independent_vars]).fit()\n",
    "\n",
    "# Print the summary results\n",
    "print(guns_model.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "66de8bf3-e8ec-4ea4-b998-194a9c518335",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                  crime   R-squared:                       0.219\n",
      "Model:                            OLS   Adj. R-squared:                  0.219\n",
      "Method:                 Least Squares   F-statistic:                     947.0\n",
      "Date:                Tue, 23 Apr 2024   Prob (F-statistic):               0.00\n",
      "Time:                        08:06:07   Log-Likelihood:                 1349.7\n",
      "No. Observations:               10144   AIC:                            -2691.\n",
      "Df Residuals:                   10140   BIC:                            -2663.\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "period         0.1192      0.006     19.443      0.000       0.107       0.131\n",
      "Israel         0.2291      0.006     38.386      0.000       0.217       0.241\n",
      "Ukraine        0.0821      0.006     14.673      0.000       0.071       0.093\n",
      "intercept      0.5575      0.005    103.664      0.000       0.547       0.568\n",
      "==============================================================================\n",
      "Omnibus:                      860.671   Durbin-Watson:                   1.939\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1092.155\n",
      "Skew:                          -0.787   Prob(JB):                    6.94e-238\n",
      "Kurtosis:                       3.326   Cond. No.                         5.35\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# Assuming df is your DataFrame\n",
    "# Add a constant term for the intercept\n",
    "df['intercept'] = 1\n",
    "\n",
    "# Specify the model\n",
    "independent_vars = ['period', 'Israel', 'Ukraine', 'intercept']\n",
    "dependent_var = 'crime'\n",
    "\n",
    "# Fit the model\n",
    "crime_model = sm.OLS(df[dependent_var], df[independent_vars]).fit()\n",
    "\n",
    "# Print the summary results\n",
    "print(crime_model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ba6bbd5",
   "metadata": {},
   "source": [
    "## Figures for paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "c65802ec-b6ad-47ec-97fa-7b1ed2a2023c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAB9CAYAAACvSWTxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAXqklEQVR4nO3de1BU5/nA8e+uLMrF9QIqirIaAgiIilHAC0ZNx9pEJo0dUybVJhGxJk2MWm20o7HTmjFWRU069TJgp6lJTdskJkZCajTWajUOiVP5ISyCCgYvCEhWWGAX9/z+IKxuUNxdlut5PjPMsO++7znveXg5z573XFajKIqCEEII1dJ2dAeEEEJ0LEkEQgihcpIIhBBC5SQRCCGEykkiEEIIlZNEIIQQKieJQAghVE4SgRBCqJxXR3egMzhz5gyKoqDT6Tq6K0II4RFWqxWNRkNsbOwD68oRAaAoCu7cYK0oChaLxa223ZHE4w6JhSOJxx3tFQtX9mtyRAD2I4GYmBiX2pnNZvLy8nj44Yfx9fVti651KRKPOyQWjiQed7RXLHJycpyuK0cEQgihcnJE4K6aGnwGDmSszUZ9SQmo/FOOEKLrkkTQChqzmR4d3QkhhGglmRoSQgiVk0QghBAqJ4lACCFUThKBEEKonCQCIYRQOblqyF1aLbcTEzHX1NBDK/lUCNF1SSJwl48P9VlZFOTlEenj09G9EUIIt8lHWSGEUDk5IhBCqF5aWhomkwm9Xs/y5cvdXo6iKBhLbvLl/12jutaKv4+O+FFBRIT0Q6PRdGjfWiKJwF01NfgMH87ohgasBQXyiAkhurC0tDRKS0sJDg52e2dbfM3Etn1nOF9cQUXpOeprbtLTrx//CI4izBDA0uRYDEH6Dunbg0giaAVNeTk6wNrRHRFCdKjiayZe/eNxCs8eI/eLdOqqK+zv9fIP4Nr0hVwtr2HjS1MYoO98u12PniO4desWv/3tb4mPj+eRRx7h9ddf56OPPiI6OhqLxcLBgweJjIyktrbWoV1KSgrLli2zv3799dd58sknOXnyJMnJyYwZM4ZZs2aRnZ3t0O7kyZPMnz+f+Ph4YmNjSUpKYv/+/Z7cJCGEaJGiKGzbd4bCs8f46sBGhyQAUFddwVcHNlJ49hjb953plN/J4LHUZLVaSU1NpaysjJUrVzJw4EB2797N4cOHGTFiBN7e3uTl5RESEoLP966yMRqNzJs3z+F1VVUVmzZtYsGCBfTu3ZvNmzezZs0asrKyAMjOziY1NZXnnnuOF154AavVSk5ODl5e7m+SzWZzpbI9i9psNtfadlNNMZB4SCy+ryvFw9X+GUsquVRaSe7R9BbrnTuawbCIeAovf2tfj6vraqvYeSwR/PnPf+bcuXNkZmYydOhQAEJDQ5kxYwZJSUkA5OXlERER4dCuoqKCGzduOJQXFBTQt29f3nnnHXvSuHHjBmvXrrXXOXDgAPHx8axYscJe9uijj7rdf4vFwrlz55yurzGbif7u9+LiYpQbN9xed3fSs2dPLl261NHd6BQkFo46czys1sYJ3qtXrzJ48GCX21ssFupuVbVYp/ZWOUfTFzDzbW+Xll1eXm7voyv7KIvFgre3c+vySCJQFIW//vWvzJs3z54EAIKDg/Hy8iIyMhJo/KSfnJzs0DY/Px/AngiuX7/OzZs3Wb58ucORw82bN+nXr5/9tY+PD6dPn2b37t0kJSW59ccTQoi72Ww2ysrK2mz5VVVVbbbs1vBIIigsLKSsrIwZM2Y4lJeXl9PQ0EBkZKT9k//IkSMd6uTn56PX6xkyZAjQmCwAJk2a1KxeWFiY/fXLL7+Moijs2bOHtLQ0xo0bx4oVKxg3bpxb2+Dt7U1UVJTzDWpq7L8aDAZ8Bwxwa73didlsxmg0EhERIV9HKLFw0Nnj0fR1tVqt1uUPlea6BqprarHWmR68nl56/P16oesBXl5eTl1SevXqVWw2GzqdzqV9VG5urtN1PZIImjJoQECAQ/np06cBGDlypP2Tf3h4eLM6d08LNSWGu48soDFBTJ482f7az8+P1atXs2rVKrKzs/nd737HkiVLOH78uNvboXXlURFeXtweN466ujq0Xl6ute2mmmKg1WpVHw+JhaOuEo/BgwfzzTffuNQmv7iSFduOcjh9UbMTxXfr1TuQx1J28friiTRUXyEyMtKppDh06FBKS0sBF/dRLvDIUpumbIqLi+1lFouFnTt3MmjQIPr3709FRWOABtz1ybmoqIj//ve/DonAaDQ2O2qwWCxcvHjRPsV0N41Gw4QJE5g5cyYNDQ2e2Bzn+PhQ/5//kP/22yCPmBBCtSJC+hFmCCB6+sIW60VPSyHcEMDDQ12/l6CteeSIICwsjODgYDZs2IDVasXLy4uMjAyuX79ObGwsAMOHDwfggw8+4IknniAnJ4etW7disVianShOSEhwWP758+dpaGiwJ4i1a9ei1WqJj48nICCAnJwcMjIyWLiw5T+EEEJ4mkajYWlyLFfLG6eLm91H0DuQ6GkpPDx6Kq8kx7p9h3Fb8kgi0Ol0bN++nXXr1vGrX/0Kg8FAamoqO3bsYPz48QDExMSwePFitm/fztatW5k0aRJLlizhF7/4hT0RNH3yf/755x2Wn5+fj06nIzQ0FGi8GungwYNkZmbS0NDAiBEjWLduHU899ZQnNkcIIVxiCNKz8aUpbNvnR1BonMOdxQHBUYQbAnjluzuLzWZzR3e3GY3SRnc3ZGVlsXLlSg4dOkRQUFBbrMJjcnJygMZk5TSzGVtkJFarldtnz+IbGNhGves6zGYzeXl5Ts99dmcSC0edPR6efNZQQclNTt31rKGEUUGE3/WsIVdj4W7fXNmveeSIIDs7m9OnTxMTE4OiKJw8eZJ33nmHZcuWdfok4DZFQVtSQk/A3AnvFBRCOM9Tz/DRaDREGPoTYejvkeWB5/rWEo8kgtraWg4ePMiuXbuAxnsC3njjDR5//HFPLF4IIUQb8kgiSExMJDEx0ROLEkII0c467wW9Qggh2oUkAiGEUDlJBEIIoXKd7xsSugqNBltkJPX19dAJbxARQghnSSJwl68vddnZjdcDd8LrooUQwlkyNSSEEConiUAIIVROpobcZTbTa/x4ourrUb78EmR6SAjRRUkicJeioM3Lwwd5xIQQomuTqSEhhFA5SQRCCKFykgiEEELlJBEIIYTKSSIQQgiVk6uG3KXRYAsJwWq1yiMmhBBdmiQCd/n6UpeXJ4+YEEJ0eTI1JIQQKieJQAghVE6mhtxVW0vPxERG1tXBsWPyiAkhRJclicBdNhs9vv4aP8Bss3V0b4QQwm0yNSSEEConiUAIIVROEoEQQqhcu54jWLhwIXq9nrS0tPZcrRCdVlpaGiaTCb1ez/Lly51upygKxpKbfPl/16iuteLvoyN+VBARIf3QuHCDo7vrF91LuyaC5cuXExAQ0J6rFKJTS0tLo7S0lODgYKd3xMXXTGzbd4bzxRVUlJ6jvuYmPf368Y/gKMIMASxNjsUQpG+z9Yvup10SgdVqRafTERUV1R6razdKYCANDQ0d3Q2hIsXXTLz6x+MUnj1G7hfp1FVX2N/r5R/AtekLuVpew8aXpjidDIRw6xyBzWZj7969PPnkk4wePZqEhAQWL15MVVUVGzduZPbs2Rw5coQ5c+YwatQoPvvsM7744gsiIiKorKwEoKamhpEjR/LPf/6T9evXk5CQQFxcHLt37wYgMzOTpKQkxo4dy/PPP29v16S8vJx169aRmJhIbGws8+bNo6CgoJXhcIGfH7XFxZz9/HPw82u/9QrVUhSFbfvOUHj2GF8d2OiQBADqqiv46sBGCs8eY/u+MyjyzXnCSS4fESiKwsqVKzly5AgpKSmMHTuWqqoqPvzwQ3x9fTEajVRWVrJlyxYWLVpEYGAg4eHhvP/++wwcOJD+/fsDYDQaURSFnTt3kpSURFpaGu+//z5btmyhrKyM8+fPs3TpUr799lvWrVvHnj17WLFiBQAVFRU8/fTT+Pn58etf/xq9Xk9GRgYLFy4kKysLXzdv7rK5eD9AU32bzeZy2+5I4nGHO7F4UD1jSSWXSivJPZreYr1zRzMYFhGPsbiC8JD+znXYifW3hoyNOzpjLFxOBO+++y7/+te/+Nvf/saoUaPs5bNnzwYgPz8fnU7Hvn376N27t/39/Px8Ro4c6fAa4Je//CVPPfUUAAaDgU8++YTc3Fz27t1Ljx49APj0008pKiqyt123bh09evTg3Xffta9jzJgxTJkyhePHjzNz5kxXNwuLxcK5c+dcbtezZ08uXbrkcrvuSuJxhzOxsFqtAFy9epXBgwc/cJkWi4W6W1Ut1qm9Vc7R9AU8+rb3A5dXXl5u74c7498VMjbuaI9YWCwWvL0fPAbAxUSgKAq7du3ipz/9qUMSaFJeXk5FRQWvvfaaQxKAxh3/D37wA4fXQ4YMsScBgNraWgAWL15sTwLQOI1kMBgAuHTpEocOHWL9+vX4+PjY5+j9/f0JDAzkypUrrmyS2zR1dQx/4YXGPu3YgdKrV7usV3RPNpuNsrIyjy2vqqrKY8sS3Z9LiaCwsJDr16/f9xN306f8qVOnOpTX1dVRUlLicERQUFBAQkJCs/ZarZa4uDiH8vPnz9vXeeLECQDWrFnDmjVrmvXBz835em9vb9dOZtfUoM3OBsAwbBi+Awa4td7uxGw2YzQaiYiIcHt6rrtwNhY6nQ4ArVb7wCMCc10D1TW1WOtMD1y/rpcefz8ffHu1/C9+9epVbDZbm1/MIWPjjvaKRW5urtN1XUoEN27cAGDAfXZ6+fn59OnTh2HDhjmUFxQUcPv2bXsiUBQFo9HIj370o2bthw8fjo+Pj72stLQUk8lEREQEAGVlZQwePJi33nrrnn1oOnJwh1brwrnzu+pqtVrX2nZTTTGQeLgei8GDB/PNN9+0WCe/uJIV245yOH1RsxPFd+vVO5DHUnaxZek0IgwtnyMYOnQopaWlDn1uCzI27uiMsXApEQQGBgKNRwYjRoxo9n5+fr59h303o9GIt7c3w4cPB6CkpASz2UxkZGSz9t8vMxqNAPbygQMHYjKZCAsLo5dMxwgViQjpR5ih8RLRrw5svG+96GkphBsCCA/p1469E12ZS4kgNDSUyMhI1q9fj8lkIjg4mMuXL/O///2P9evXYzQamThxYrN2RqORsLAwvLy87K8Bh6kiaEwEP//5z5u1DQoKom/fvgA89thjbN26ldTUVObPn0/fvn0pLy8nOzubqVOnMm3aNFc2SYguQ6PRsDQ5lqvlNQDN7yPoHUj0tBQeHj2VV5JjXbrDWKibS4mgR48e/OlPf+IPf/gDaWlpVFdXM3ToUJ555hksFgsXL15kwYIFzdp9/0ghPz+f4OBg9Po7N7xUVlZy48aNex4R3F0WFBTEnj172LZtG2vXrqW+vp5BgwYRFxdHdHS0K5sjRJdjCNKz8aUpbNvnR1BonMOdxQHBUYQbAnjFhTuLhQDQKHLXCTk5OQDExMQ436imBvz9ATCXlcnJYhpPguXl5REZGSknBJ2MRWueNVRQcpNTdz1rKGFUEOGd9FlDMjbuaK9YuLJfky+maQXF17fT3BAiuiZ3d74ajYYIQ/8Hngxuq/WL7kUSgbv8/Ki9caMxs8sjJoQQXVjnuHZJCCFEh5FEIIQQKidTQ+6qq6PnnDk8XF0NH38MKj8BJoTouiQRuOv2bXp89hl9APPt2x3dGyGEcJtMDQkhhMpJIhBCCJWTRCCEEConiUAIIVROEoEQQqicPGsI+Prrr1EUxemvdQNAUeDixcZfDQY0d32jmlopioLVakWn06n+yZcSC0cSjzvaKxYWiwWNRsO4ceMeWFcuHwX3/hgaDTz0UOOvHu5PV6XRaFxLpt2YxMKRxOOO9oqFRqNxet8mRwRCCKFyco5ACCFUThKBEEKonCQCIYRQOUkEQgihcpIIhBBC5SQRCCGEykkiEEIIlZNEIIQQKieJQAghVE4SgRBCqJwkAiGEUDlJBMDFixdJSUlh7NixTJw4kfXr11NXV+dU2w8//JBZs2YRExPD7Nmz+fTTT5vVsVqtbNmyhSlTpjBmzBjmz59Pfn6+pzfDY9o6HhEREc1+Jk+e7OnN8Ah3Y5GZmcnLL79MYmIiERERZGRk3LOeWsaGs/Ho7mOjurqat956i7lz5zJ+/HgSEhJISUkhNze3Wd32HBuqf/qoyWTi2WefZciQIbz55ptUVlayYcMGqqqq2Lx5c4tts7KyWLVqFYsWLWLy5Ml8/vnnLFu2jN69ezNlyhR7vQ0bNrB//35WrVpFcHAw6enpPPfccxw4cIABAwa09Sa6pD3iATB//nxmz55tf63T6dpke1qjtbG4fPky06dP57333rtvPTWNDWfiAd17bFy5coX33nuPn/zkJyxZsoSGhgbefvttkpOT2bdvH9HR0fa67To2FJXbtWuXMmbMGKWiosJe9vHHHyvh4eFKYWFhi21nzZqlLFmyxKFswYIFyty5c+2vr127pkRGRip79+61l926dUuJi4tTNm3a5KGt8Jy2joeiKEp4eLiSnp7uuU63kdbE4vbt2/bf77e9ahobzsTjQe91Ju7GoqamRjGbzQ5ldXV1yuTJk5VVq1bZy9p7bKh+aujYsWNMnDiR/v3728t++MMf4u3tzb///e/7trt8+TIXLlxw+OQCMHv2bM6ePUtlZSUAx48f5/bt2zzxxBP2Ov7+/syYMaPF5XeUto5HV+JuLAC02gf/a6llbIBz8ehK3I2Fr68vPj4+DmU9e/YkNDSUsrIye1l7j43u9ddxQ1FREaGhoQ5l3t7ehISEUFRUdN92Fy5cAOCh776cpkloaCiKotjfLyoqIjAwkL59+zard/HiRWw2mwe2wnPaOh5Ndu/eTXR0NOPHj2fp0qVcuXLFQ1vgOe7GwpXlq2FsuEptY8NsNpOXl+fwv9PeY0POEZhM6PX6ZuV6vZ5vv/32vu2a3vt+2z59+ji8bzKZ6N27d7P2ffr0wWq1Yjab8ff3d7v/ntbW8QD48Y9/zLRp0wgMDKSgoIAdO3bwzDPP8NFHH9nrdwbuxsKV5athbLhCjWNj27Zt1NbWMm/ePIflt+fYUH0iuB9FUZz6mrfv11G++8K3u8vvtRyli30xnCfjsXHjRvvvEyZM4JFHHmHOnDn8/e9/JzU11UM9bjvOxsIZahobzlDb2Dhw4AB/+ctfeO211zAYDA7vtefYUP3UkF6vx2QyNSu/devWPTN+k3t90gXsy2pqe7/lm0wmdDodvr6+bve9LbR1PO5l5MiRjBgx4p6X0HUkd2PR2uV3t7HRGt15bJw4cYLVq1eTkpLCz372M6eW31ZjQ/WJIDQ0tNmcnsVioaSkpNkc4N2a5vO+P/ddVFSERqOxvx8aGkpFRQVVVVXN6o0YMaLTnURr63jcT2f8FOxuLFxZvhrGRmt1x7Fx9uxZXnrpJWbNmsXKlSvvufz2HBuda6R1gKlTp3Lq1Clu3rxpLzt06BAWi4VHH330vu2GDRvGQw89RGZmpkP5J598wujRo+1XE0yZMgWtVutwY1VNTQ1Hjhxpcfkdpa3jcS95eXlcunSJmJiY1m+AB7kbC2epZWy0RnccG0VFRaSmpjJu3Dg2bNhwzymg9h4bqj9HkJyczN69e3nxxRd58cUXqaio4I033iApKckhs//mN79h//79nDt3zl62ZMkSli1bRkhICJMmTeLw4cOcOHGC9PR0e51BgwaRnJzM5s2b8fLyYsiQIezZsweAZ599tv021EltHY+MjAwuX75MXFwc/fv35/z58+zcuZOgoCDmzp3brtv6IK2JRWFhIYWFhfbXBQUFZGVl4ePjY/9HVtPYcCYeahgbFRUVpKSkoNPpWLhwocOUl7e3N1FRUUAHjA2P35nQBV24cEFZsGCBMmbMGCU+Pl75/e9/r9TW1jrUefXVV5Xw8PBmbT/44ANl5syZSnR0tPL4448rmZmZzerU19crmzZtUiZNmqTExMQo8+bNU/Ly8tpse1qrLeNx+PBh5emnn1YmTJigREVFKZMnT1ZWr16tXL9+vU23yV3uxuLNN99UwsPDm/1Mnz7doZ5axoYz8VDD2Dh16tQ949DRY0OjKJ1wAk4IIUS7Uf05AiGEUDtJBEIIoXKSCIQQQuUkEQghhMpJIhBCCJWTRCCEEConiUAIIVROEoEQQqicJAIhhFA5SQRCCKFykgiEEELlJBEIIYTK/T9UmOvy/Ast1QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x100 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get forest plot showing the confidence interval for the coefficient of the period variables for guns and crime\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "# Create a DataFrame with four columns: Measure (the variable names), Estimate (the coefficients), and CI_lower, CI_upper (the confidence intervals)\n",
    "forest_df = pd.DataFrame({'Measure': ['crime', 'guns'],\n",
    "                            'Estimate': [crime_model.params['period'], guns_model.params['period']],\n",
    "                            'CI_lower': [crime_model.conf_int().loc['period'][0], guns_model.conf_int().loc['period'][0]],\n",
    "                            'CI_upper': [crime_model.conf_int().loc['period'][1], guns_model.conf_int().loc['period'][1]]})\n",
    "\n",
    "# Plot the forest plot. Include a dashed red vertical line at 0 to indicate the null hypothesis.\n",
    "plt.figure(figsize=(4, 1))\n",
    "sns.pointplot(data=forest_df, x='Estimate', y='Measure', linestyle='none')\n",
    "plt.axhline(y=0, color='lightgrey', linestyle='-')  # Add horizontal grid line\n",
    "plt.axhline(y=1, color='lightgrey', linestyle='-')  # Add horizontal grid line\n",
    "plt.errorbar(x=forest_df['Estimate'], y=forest_df['Measure'], xerr=[forest_df['Estimate'] - forest_df['CI_lower'], forest_df['CI_upper'] - forest_df['Estimate']], fmt='o', color='black', capsize=4, linewidth=2, capthick=2)\n",
    "plt.axvline(x=0, color='red', linestyle='--')\n",
    "plt.ylabel('')  # Remove y-axis label\n",
    "plt.xlabel('')  # Remove x-axis label\n",
    "\n",
    "# Adjust y-axis limits\n",
    "plt.ylim(-0.5, 1.5)\n",
    "\n",
    "# Make y-axis tick labels italic\n",
    "plt.yticks(fontstyle='italic')\n",
    "\n",
    "# Save the plot as a PNG file with high dpi\n",
    "plt.savefig('fig_3.png', dpi=300, bbox_inches='tight')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "catching_trolls",
   "language": "python",
   "name": "catching_trolls"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
